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ABSTRACT 

The Eulerian cosmological fluid equations are used to study the nonhnear 
mode couphng of density fluctuations. We evaluate the second-order power 
spectrum including all four-point contributions. In the weakly nonlinear regime 
we find that the dominant nonlinear contribution for realistic cosmological 
spectra is made by the coupling of long-wave modes and is well estimated by 
second order perturbation theory. For a linear spectrum like that of the cold 
dark matter model, second order effects cause a significant enhancement of the 
high k part of the spectrum and a slight suppression at low k near the peak of 
the spectrum. Our perturbative results agree well in the quasilinear regime with 
the nonlinear spectrum from high-resolution N-body simulations. 

We find that due to the long-wave mode coupling, characteristic nonlinear 
masses grow less slowly in time (i.e., are larger at higher redshifts) than would 
be estimated using the linear power spectrum. For the cold dark matter model 
at {1 + z) = (20,10,5,2) the nonlinear mass is about (180,8,2.5,1.6) times 
(respectively) larger than a linear extrapolation would indicate, if the condition 
rms Sp/p = 1 is used to define the nonlinear scale. At high redshift the 
Press-Schechter mass distribution significantly underestimates the abundance of 
high-mass objects for the cold dark matter model. Although the quantitative 
results depend on the definition of the nonlinear scale, these basic consequences 
hold for any initial spectrum whose post-recombination spectral index n 
decreases sufficiently rapidly with increasing k, a feature which arises quite 
generally during the transition from a radiation- to matter-dominated universe. 



Subject headings: cosmology: theory — large-scale structure of universe — 
galaxies: clustering — galaxies: formation 
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1. Introduction 

There exists a standard paradigm for the formation of cosmic structure: gravitational 
instabihty in an expanding universe. According to this paradigm, dark matter density 
fluctuations 5{x ) = 5p{x ) / p created in the early universe lay dormant until the universe 
became matter-dominated at a redshift lO'^O/i^ (where VL is the 

present density parameter for nonrelativistic matter and the present Hubble parameter is 
Hq = 100 h kms^^ Mpc^^). After this time, the density fluctuations increased in amplitude 
as predicted by the well-known results of linear perturbation theory (e.g., Peebles 1980; 
Efstathiou 1990; Bertschinger 1992), until the fluctuations became nonlinear on some length 
scale. Bound condensations of this scale then collapsed and virialized, forming the first 
generation of objects (Gunn & Gott 1972; Press & Schechter 1974). Structure formation 
then proceeded hierarchically as density fluctuations became nonlinear on successively 
larger scales. 

At early times density fluctuations were small on the length scales of present day 
large-scale structure. Therefore, after the universe became matter dominated fluctuations 
on scales much larger than the scale of collapsed objects can be studied under the 
approximation of a pressureless, irrotational fluid evolving under the action of Newtonian 
gravity. A pcrturbative analysis of the fluid equations in Fourier space can then be used 
to study the effects of mode coupling between scales that arc weakly nonlinear. This is 
the approach we shall follow in this paper. Nonlinear analyses in real and Fourier space 
are somewhat complementary in that real space analyses are best suited to studying the 
effect of nonlinearities on the collapse and shapes of individual objects (Bertschinger & 
Jain 1993), whereas Fourier space studies provide estimates of how different parts of the 
initial spectrum couple and influence the evolution of statistical quantities like the power 
spectrum. In principle of course, the two approaches are equivalent and should give the 
same information. For pcrturbative analyses in real space see, e.g., Peebles (1980), Fry 
(1984), Hoffman (1987), Zaroubi & Hoffman (1993), and references therein. 

AlthoTigh density fluctuations of different wavelengths evolve independently in linear 
perturbation theory, higher order calculations provide an estimate of some nonlinear 
effects. Preliminary second order analyses have led to the conventional view that in models 
with decreasing amounts of power on larger scales long-wavclcngth fluctuations have no 
signiflcant effect on the gravitational instability occuring on small scales. On the other 
hand, it is known that under some circumstances small-scale, nonlinear waves can transfer 
significant amounts of power to long- wavelength, linear waves. If the initial spectrum is 
steeper than A;^ at small k (comoving wavenumber), then small-scale, nonlinear waves can 
transfer power to long wavelength linear waves so as to produce a A;^ tail in the spectrum. 
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(Zel'dovich 1965; Peebles 1980, Section 28; Vishniac 1983; Shandarin & Melott 1990). 

The question of whether power can be transfered from large to small scales was 
examined by Juszkiewicz (1981), Vishniac (1983), Juszkiewicz, Sonoda & Barrow (1984), 
and more recently by Coles (1990), Suto and Sasaki (1991) and Makino, Sasaki and Suto 
(1992). Their analyses involved writing down integral expressions for the second order 
contribution to the power spectrum, examining their limiting forms and evaluating them for 
some forms of the initial spectrum. Juszkiewicz et al. (1984) examined the autocorrelation 
function and found that the clustering length decreases due to power transfer from large to 
small scales for the initial spectrum P{k) oc A;^. However, for the cold dark matter (hereafter 
CDM) spectrum Coles (1990) found the opposite effect, though it is not significant unless 
(78 is taken larger than 1. Makino et al. (1992) have analytically obtained the second 
order contributions for power law spectra, and estimated the contribution for the CDM 
spectrum by approximating it as two power laws. Bond & Couchman (1988) have compared 
the second order CDM power spectrum to the Zel'dovich approximation evaluated at 
the same order. Some issues of mode coupling have recently been investigated through 
N-body simulations in 2-dimensions (see e.g., Beacom et al. 1991; Ryden & Gramann 1991; 
Gramann 1992). 

We have used the formalism developed in some of the perturbative studies cited 
above, and especially by Goroff et al. (1986), to calculate second order contributions to 
the power spectrum (i.e., up to fourth order in the initial density) for the standard CDM 
spectrum. Second order perturbation theory has a restricted regime of validity, because 
once the density fluctuations become sufficiently large the perturbative expansion breaks 
down. For this reason N-body simulations have been used more extensively to study the 
fully nonlinear evolution of density fluctuations. However, perturbation theory is very 
well suited to address some speciflc aspects of nonlinear evolution and to provide a better 
understanding of the physical processes involved. Being less costly and time-consuming 
than N-body simulations, it lends itself easily to the study of different models. Perturbation 
theory should be considered a complementary technique to N-body simulations, for while 
its validity is limited, it does not suffer from the resolution limits that can affect the latter. 
Hence by comparing the two techniques their domains of validity can be tested and their 
drawbacks can be better understood. In this paper we shall make such comparisons for the 
CDM spectrum. 

The most powerful use of perturbative calculations is in the study of weakly nonlinear 
evolution out to very high redshifts, spanning decades of comoving length scales in 
the spectrum. Since the formulation of the perturbative expansion allows for the time 
evolution of the spectrum to be obtained straightforwardly, we obtain the scaling in time of 
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characteristic nonlinear mass scales ranging from the nonlinear scale today, about 10^ Mq, 
to about 10^ Mq, the smallest baryonic mass scale likely to have gone nonlinear after the 
universe became matter dominated. Such an analysis cannot be done by existing N-body 
simulations as the dynamic range required to cover the full range of scales with adequate 
spectral resolution exceeds that of the current state-of-the-art. 

There are two principal limitations to our analytic treatment: the first arises from 
the general problem that the perturbative expansion breaks down when nonlinear effects 
become sufficiently strong. This drawback is particularly severe in our case because the 
regime of validity is difficult to estimate. It is reasonable to expect that second order 
perturbation theory ceases to be valid when the rms 5p/ p ^ 1, but one cannot be more 
precise without explicitly calculating higher order contributions. 

The second kind of limitation arises from the simplifying assumptions that pressure 
and vorticity are negligible. On small enough scales nonlinear evolution causes the 
intersection of particle orbits and thus generates pressure and vorticity. Through these 
effects virialization on small-scales can alter the growth of fluctuations on larger scales. It 
is plausible that the scales in the weakly nonlinear regime are large enough that this effect 
is not significant. This belief is supported by heuristic arguments as well as recent studies 
of N-body simulations (Little, Weinberg & Park 1991; Evrard & Crone 1992 and references 
therein). We conclude that the first kind of limitation, namely the neglect of higher order 
contributions, or worse still, the complete breakdown of the perturbative expansion, is likely 
to be more severe for our results. We shall address this where appropriate and accordingly 
attempt to draw conservative conclusions supported by our own N-body simulations. 

The formalism for the perturbative calculation is described in Section 2. We describe 
the numerical results for CDM in Section 3.1 and compare them to N-body simulations 
in Section 3.2. The scaling of the nonlinear scale as a function of redshift is presented in 
Section 3.3. The distribution of nonlinear masses is examined in Section 3.4 We discuss 
cosmological implications of the results in Section 4. 



2. Perturbation Theory 

In this section we describe the formalism for perturbative solutions of the cosmological 
fluid equations in Fourier space. Our approach is similar to that of Goroff et al. (1986). 
The formal perturbative solutions are then used to write down the explicit form of the 
second order contribution to the power spectrum. 
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2.1. General Formalism 



We suppose for simplicity that the matter distribution after recombination may be 
approximated as a pressureless fluid with no vorticity. We further assume that peculiar 
velocities are nonrelativistic and that the wavelengths of interest are much smaller than 
the Hubble distance cH~^ so that a nonrelativistic Newtonian treatment is valid. Using 
comoving coordinates x and conformal time dr = dt/a{t), where a{t) is the expansion scale 
factor, the nonrelativistic cosmological fluid equations are 

|^ + V-[(l + 5)t?]=0, (la) 

dv ^ a ^ -± , , 

— +[v-V]v = — lb 

= A-KGa^pS , (Ic) 

where d = da /dr. Note that v = dx/dr is the proper peculiar velocity, which we take to be 
a potential field so that v is fully specified by its divergence: 

e = W-v . (Id) 

We assume an Einstein-de Sitter {Q = 1) universe, with a oc 

^2/3 

OC r . We will also assume 
that the initial (linear) density fluctuation field is a gaussian random field. 

To quantify the amplitude of fluctuations of various scales it is preferable to work with 
the Fourier transform of the density fluctuation field, which we define as 

d^x 



S{k,T) = y (2) 

and similarly for §(k,T). The power spectrum (power spectral density) of S{x,t) is defined 
by the ensemble average two-point function, 

{6{h,r)6{k2,T)) = P{k,,T)6T){h + k2) , (3) 

where is the Dirac delta function, required for a spatially homogeneous random density 
field. For a homogeneous and isotropic random field the power spectrum depends only on 
the magnitude of the wavevector. The contribution to the variance of 6{x,t) from waves in 
the wavevector volume element d^k is P{k,T)d^k. 



Fourier transforming equations (1) gives 

^ + e=-Jd'kiJd^k2ST,{h ..^. ^2 



— + e = - J d^kij d^k2Sji{ki + k2 - k) -j-^e{ki,T) S{k2,T) , (4a) 
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^ + l^ + ^^ = -J d'kij d'h Snih + h-k) Hh: r) ^(^2, r) . (4b) 

In equations (4) the nonlinear terms constitute the right-hand side and illustrate that 
the nonlinear evolution of the fields S and ^ at a given wavevector k is determined by the 

— * 

mode coupling of the fields at all pairs of wavevectors whose sum is k, as required by spatial 
homogeneity. This makes it impossible to obtain exact solutions to the equations, so that 
the only general analytical technique for self-consistently evaluating the nonlinear terms is 
to make a perturbative expansion in 6 and 9. The formalism for such an expansion has 
been systematically developed by Goroff et al. (1986) and recently extended by Makino et 
al. (1992). Following these authors we write the solution to equations (4) as a perturbation 
series, 

00 00 
5{k, t) = Y: «"(t) 5n{k) , §{k, t)^Y1 «(^)«""'(^) Uk) (5) 

n=l ri=l 

It is easy to verify that for n = 1 the time dependent part of the solution correctly gives 
the linear growing modes 61 oc a(r) and 9i (x a and that the time-dependence is consistent 
with equations (4) for all n. To obtain formal solutions for the k dependence at all orders 
we proceed as follows. 

Substituting equation (5) into equations (4) yields, for n > 1, 

nSn{k) + 9n{k) = A^{k) , SS^ik) + (1 + 2n)9n{k) = B^{k) , (6) 

where 

A^k ) = - j d^ki J d'k2 5d(A;i + k2-k ) ^ 9m{kl) 6n-m{k2) , (7a) 



'^l m=l 

d^k^ / d^'k^ 5j,{ki + k2-k) \' ^ ^„(A;i) 9n-m{k2) . (7b) 

•' 1^11^2 rn=l 

Solving equations (6) for 5n and 9n gives, for n > 1, 

^ {l + 2n)An{k)-Bn{k) _ -3A„(fc)+n5„(^) 
(2n + 3)(n-l) ' ^^^^^ ' (2n + 3)(n-l) " 

Equations (7) and (8) give recursion relations for 5n{k) and 9n{k), with starting values 
6i{k) and 9i = —Si. The general solution may be written 

5n{k ) = j d^qi ■■■ j d^Qn Sniqi H \-qn-k) -F„(^i, . . . , g„) 5i{qi) • • • 5i{qn) , (9a) 

9n{k ) = - / ^^^^1 • • • / d^qn 5D(gi H \-qn-k) Gn{qi, ...,qn) Si{qi) ■ ■ ■ Si{qn) ■ (9b) 
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Prom equations (7)-(9) we obtain recursion relations for F„ and G„: 
Fn(qi,...,qn) 



m 



[ (2n + 3)(n- 1) 



1 + '^ri)^jP-Fn-m{qm+l: ■ ■ ■ : qn) 



Gn—m{,qm+li ■ ■ ■ i qn) 



(10a) 



Gn{qi,---,qn) 



Gjji^qij • • • ? Qm) 

m= 



+n 



\ (2n + 3)(n- 1) 
k\ki ■ h 



k-ki 



^n— m(?m+l 



Gn—m{qm+lj • • • ) 5n) 



,?n) 



(10b) 



where = gi + • • • + q^, k2 = qm+i + + ^ = + ^2 and Fi = Gi = 1. Equations 
(10) are equivalent to equations (6) and (Al) of Goroff et al. (1986), with F„ = P„ and 
Gn — (3/2)(5n in their notation. 



2.2. Power Spectrum at Second Order 



To calculate the power spectrum we shall prefer to use symmetrized forms of F„ and 
On, denoted F^^^ and Cl^^ and obtained by summing the n! permutations of F„ and G„ 
over their n arguments and dividing by n\. Since the arguments are dummy variables of 
integration the symmetrized functions can be used in equations (9) without changing the 
result. The symmetrized second-order solutions of equations (10) are given by 



Fi'\kuk2) 



5 2 {kl ■ k2f , {kl ■k2) ( I 1 



^1^2 



+ 



Ms),r 3 A {ki.k2f {h-h) 

G, iki,k2)-- + -^^ + ^^ 



:iia) 



(lib) 



Note that ^2*"* and 0^2"^ have first-order poles as ki 



or k2 for fixed k: 



G2^ ~ (1/2) cos'd {ki/k2 + k2/ki) where d is the angle between ki and k2- The 
expression for Fg"^ will also be required, but since it is very long we shall wait to write a 
simplified form below. 
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The recursion relations in equations (10) may be used to compute the power spectrum 
at any order in perturbation theory. Substituting equation (5) into equation (3), we have 



{5{k,r)5{k\r)) 
a'{r){Si{k) Si{k')) + a\r) 



{Si{k)S^{k')) + {S2{k)S2{k')) 



-iWSiik')) 



(12) 



Equation (12) exphcitly shows all the terms contributing to the power spectrum at fourth 
order in the initial density field 5i (or second order in the initial spectrum) , as the nth order 
field Sn{k) involves n powers of Si{k). With the definition 



{Sm{k) 5n^rn{k!)) = P„,„_„(A;) Sj){k + k') (13) 

the power spectrum up to second order (i.e., fourth order in Si) is given by equation (12) as 

P{k,T) = a^{T)Pn{k) + a''{T)[P22{k) + 2Pis{k)] 

= a'{T)Pn{k) + a\r)P2{k), (14) 

where the net second order contribution P2{k) is defined as 

P2{k)^P22{k)+2Pisik). (15) 

To determine P2{k) we need to evaluate the 4-point correlations of the linear density 
field 6i{k). For a gaussian random field, all cumulants (irreducible correlation functions) 
of Si{k) vanish aside from the 2-point cumulant, which is given by equation (3) for 
m — n — m — 1. All odd moments of di{k ) vanish. Even moments are given by symmetrized 

— * 

products of the 2-point cumulants. Thus the 4-point correlation function of 5i{k) is 



(Siih) Siih) Siih) 5i{h)) = P{ki)P{ks)6j,ih + k2)6j,ih + h) 
+P{ki)P{k2)6T,{h + 4)5d(4 + h) + P{ki)P{k2)6j,{h + h)5j,{k2 + h) 



(16) 



With the results and techniques described above, we can proceed to obtain the second 
order contribution to the power spectrum. The two terms contributing at second order 
simplify to the following 3-dimensional integrals in wavevector space: 

2 



P22{k) = 2jd\ Piiiq) Pii{\k-q\) [F^'\q, k - q) 



with ^2*-* given by equation (11a), and 



2Pis{k) = 6Pn{k) J d'qPniq) F^'\q, -q,k) 



;i7) 



(18) 



9 



The numbers in front of the integrals arise from the procedure of taking expectation values 
illustrated in equation (16). We write the integrals in spherical coordinates g, 'i?, and (p: 
the magnitude, polar angle and azimuthal angle, respectively, of the wavevector q. Then 
with the external wavevector k aligned along the z-axis the integral over is trivial and 
simplifies / d^q to 27i J dqq^ J dcosi). For P13, the dependence on i) is also straightforward 
as it arises only through F-^^^ and not Pn. This allows the integral over cos'd to be done 
analytically as well, giving (Makino et al. 1992) 



Thus with a specified initial spectrum Pu{k) equations (17) and (19) give the second 
order contribution. Before evaluating these integrals for the CDM initial spectrum, we 
point out that the poles of F2 and G2 described after equations (11) give the leading order 
part of the integrand of equation (17) in (q/k) as: 



If Pu{k) ~ A;" with n < — 1 as /c — * 0, then P22 diverges. Vishniac (1983) showed that the 
leading order part of 2P13 in {q/k) is negative and exactly cancels that of P22 — this can be 
demonstrated by examining the limiting form of F^^\ In a future paper we will analyze the 
leading order behavior of perturbative integrals at higher orders and also calculate it using 
a nonperturbative approach in order to investigate whether there may exist divergences for 
some power spectra at higher orders in perturbation theory. For the purposes of the second 
order integration the cancellation of the leading order terms has no consequence other than 
requiring that each piece, P22 and P13, be integrated very accurately to get the resultant. 
This is necessary because the cancelling parts cannot be removed before performing the 
integrals as the two integrands have different forms: P22 is symmetric in q and {k — q), 
whereas P13 is not. We will return to this point in the next section. 



The results obtained in the previous section will now be used to obtain the second 
order contributions to the CDM power spectrum. We will use the standard CDM spectrum 




(19) 




(20) 



3. 



Results for CDM 
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with parameters Q — 1, Hq — 50 kms ^Mpc ^, and as — 1. For the linear spectrum at 
a = 1 we use the fitting form given by Bardeen et al. (1986): 



Puik) = AkT^{k) , A = 2.19 x 10%Ipc^ 



ln(l + 9.36/c) 
9.36fc 



1 + 15M + {64:Akf + {21.8kf + (26.8A;)' 



-1/4 



(21) 



where k is in units of Mpc^^. With this initial spectrum equations (17) and (19) can be 
used to obtain the second order contribution P2{k), which can then be used to obtain the 
net power spectrum as a function of a and k from equation (14). 



3.1. Nonlinear Power Spectrum 



As pointed out in Section 2.2 the integrals for P22 and P13 contain large contributions 
which exactly cancel each other. For the CDM spectrum these contributions are finite but 
care is still required in their numerical evaluation. Equal contributions from P22 are made 
as g* — > and q ^ k, whereas the cancelling contribution from 2P13 is made only as if 0. 
The integrand for P22 is symmetric in q and {k — q) and is positive definite. For ease of 
numerical integration, we break up the integration range for P22 as follows: 

(ff'q re j-l rk-t r-l rk+e Ak^ -e^)/2kq 

— — — 2 dq dy + dq dy + dq dy 



27r JO J-l Je J-l Jk-e J-l 

fk-kc fl fkc fl 

+ dq dy+ dq dy , (22) 

Jk+t J-l Jk-kc J(k^W-kl)/2kq 

where y = cosi), and kc is the upper limit required because at high q the spectrum has 
departed strongly from the linear spectrum causing the perturbative expansion to break 
down. Transfer of power from higher frequencies is suppressed by virialization. The first 
term on the right-hand side of equation (22) has a factor of 2 because we have used the 
symmetry between q and {k — q) in the integrand to exclude a small ball of radius e 

— * 

around q — k (where the integration becomes difficult) by restricting the limits on y in 
the third term, requiring us to double the contribution from a similar ball around q — 

— * 

to compensate. The limits on y in the last term are set to ensure that \k — q\ < kc as 
required to consistently impose the upper limit, i.e., to exclude any contribution from Pn 
in equation (17) when its argument exceeds kc- It is in principle important to scale kc with 
time to reflect the growth of the nonlinear length scale with time, because that determines 
the range of validity of the perturbative expansion. We have done so using the linear scaling 
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ke (X a 2/(3+"), although as explained below at early times the result is insensitive to the 
choice of kc- 

The results of performing the integrals in equations (17) and (19) for a large range of 
values of k are shown in Figure 1. We plot the hnear spectrum a^Pii(A;), the net spectrum 
including second order contributions given by equation (14), and the nonlinear spectrum 
computed from high-resolution N-body simulations described in Section 3.2 at four values 
of the expansion factor. The spectra have been divided by to facilitate comparison of 
the results at different times. The second order results at different values of a are obtained 
by simply multiplying Pu and P2 by different powers of a as shown in equation (14), so 
the integration of P22 and P13 needs to be done only once for a given k. The second order 
spectrum should be taken seriously only for the range of k for which a^P2{k) < a'^Pn{k), 
as we do not expect the perturbative results to be valid for higher k. The interesting range 
of k, the regime where nonlinear effects set in, moves to lower k as one looks at larger a, 
reflecting the progress of nonlinearities to larger length scales (lower k) at late times. As 
expected we find that at a given time the second order contribution is not significant for 
small k where the rms Sp/p 1. 

For small k up to just over the peak of the spectrum, the second order contribution is 
negative, causing the nonlinear spectrum to be lower than the linear one. At relatively high 
k the second order contribution enhances the growth of the spectrum. This has the effect 
of making the slope of the spectrum significantly shallower at high k than that of the linear 
spectrum. Thus, power is effectively transfered from long to short wavelengths, although 
the enhancement at short wavelengths exceeds the suppression at long wavelengths. 

The two power law model of Makino et al. (1992) gives qualitatively similar results 
to those shown in Figure 1. Bond & Couchman (1988) also computed the second order 
contributions to the CDM spectrum with a view to checking the reliability of the Zel'dovich 
approximation at the same order. They found excellent agreement, in contrast to the 
results of Grinstein & Wise (1987) who found that in comparison to perturbation theory 
the Zel'dovich approximation significantly underestimated the magnitudes of the gaussian 
filtered, connected parts of the third and fourth moment of the real space density. In 
comparison to our results. Figure 3 of Bond & Couchman shows a larger enhancement over 
the linear spectrum, and does not appear to show the suppression at relatively low k at 
all. They do not give the explicit form of the term corresponding to our P13, but state that 
it is negligible in comparison to P22. This does not agree with our results at low k and is 
probably the source of the difference in our figures. It is difficult to make a more detailed 
comparison without knowing the explicit form of their second term. 

In order to obtain a better understanding of the dynamics of the mode-coupling, we 
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have examined the relative contribution of different parts of the CDM initial spectrum to 

— * 

the second order results at a given k. Let q denote the integrated wavevector and k the 
external wavevector at which the second order contribution is calculated, as in equations 
(17) and (19). There is a two-fold ambiguity because wavevector k — q contributes at the 
same time as q. We have carefully examined different ways of associating second order 
contributions from different parts of the initial spectrum, and found that the second order 
contribution from q ^ k tends to be positive and that from q ^ k negative. Indeed we also 
find this to hold for power law spectra with —3 < n < 1, independently of the value of n, 
thus indicating that it is a general feature of second-order mode-coupling. In Figure 2 we 

— * 

associate the second-order contributions, dP2{k) / d\uq, with the smaller of q and |A; — ^|. 

There are two regimes in the CDM spectrum, divided roughly by the part where the 
logarithmic slope n[= d\n.Pii)/dhik\ falls below —1. For small k, where n ^ — 1, the 
positive second order contribution from ^ ^ A; is swamped by the negative contribution 
made by large q. The net effect is to decrease the growth of the spectrum compared to the 
linear growth. For relatively large k, where n ^ —2 the positive contribution from small q 
dominates, increasingly so as one goes to higher k. A comparison of the curves in Figure 
2 for A; = 0.1 Mpc^-*^ and k = 1 Mpc^^ shows how the relative strengths of the positive and 
negative contributions shift as one moves across the spectrum. This shift can be understood 
by observing that at higher k there is an increasing amount of power in the initial spectrum 
a,t q < k; the plot of the rms power on scale k in Figure 3 illustrates this point. The 
increased power at small q causes a larger nonlinear enhancement at higher k. Since the 
weakly nonlinear regime moves to higher k at earlier times, the enhancement at high k 
in turn leads to a stronger nonlinear growth at earlier times. We study the consequence 
of this fact in detail in Section 3.3. The dominance of the nonlinear contribution from 
long-wave modes also strengthens the consistency of the perturbative calculation, because 
the amplitude of the density fluctuations is small for these modes. As discussed in Section 
3.2, this may be responsible for the second order results being valid for a much larger range 
of scales at earher times. 

We emphasize that the transition value n ~ — 1 for the change in sign of the second 
order contribution is only approximate, because it depends on the value of k taken as 
being representative of the weakly nonlinear regime. We have examined the second order 
contribution for power law spectra Pii{k) oc A;", for a range of values of n between —3 
and 1 to verify this transition. We find that in the weakly nonlinear regime (defined by 
k ^ kni, where kni is the scale at which the rms 5p/p — 1), the second order contribution 
for n sufficiently larger than —1 is negative and that for n sufficiently smaller than —1 is 
positive. For n ~ — 1, the contribution is negative for low k and positive for high k in the 
weakly nonlinear regime. These results are consistent with the results of Makino et al. 
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(1992) who examined the second order contributions for n = 1, 0, —1, —2; they also found 
good agreement with N-body simulations. 

A possible transition in the nature of nonlinear evolution at n = — 1 has also been 
explored by studying the clustering in real space in N-body simulations by Klypin & Melott 
(1992). An examination of the origin of the term providing the dominant second order 
enhancement suggests that the advective {v ■ V) terms in the real space fluid equations 
cause the change in sign of the nonlinear contribution. This interpretation is consistent 
with the fact that for n < —1 there is an increasing amount of power in the rms velocity 
field on larger scales, and this appears to cause the nonlinear enhancement of the density 
from long-wave modes to dominate. These arguments are by no means rigorous, and merit 
further exploration. 

It is worth noting that for the deeply nonlinear regime the stable clustering hypothesis 
(Peebles 1980, Section 73) predicts that the spectrum steepens below the linear theory 
spectrum for n > —2 and rises above it for n < —2. Consistency with the second order 
results would require that at least for — 2 ^ n ;S — 1, the nonlinear spectrum first rise 
above the linear one in the weakly nonlinear regime and then fall below it in the deeply 
nonlinear regime. This is indeed seen in N-body simulations; the results of Efstathiou et al. 
(1988) show only hints of this feature owing to limited resolution, but it is clearly evident 
in simulations with higher resolution (Bertschinger & Gelb 1991; White 1993). 

The second order spectrum provides an estimate of the change in the fiuctuation 
amplitude due to nonlinear effects in the weakly nonlinear regime. The conventional 
normalization is to set the rms 5p/p on a scale of 8h~^ Mpc, denoted as, equal to 1. The rms 
value is computed from the power spectrum using a top-hat filter as described in Section 
3.3. We find that with the linear spectrum normalized in this way, second order effects 
increase ag by 10%. This is a smaller enhancement than found by Hoffman (1987) for the 
standard deviation of the density (without filtering) using the Zel'dovich approximation. 
The N-body spectrum shows an even smaller change in erg than the second order spectrum, 
although it is difficult to estimate accurately in a box of length 50 Mpc. 



3.2. Comparison with N-Body Simulations 

The N-body results shown in Figure 1 are from two different particle-particle/particle- 
mcsh simulations of the CDM model in a (100 Mpc)^ box normalized so that linear ag = 1 
at a = 1/ {1 + z) = 1. For a > 0.1 we have used the simulation with 144^ particles and 
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Plummer softening distance 65 kpc performed by Gelb & Bertschinger (1993). To obtain 

accurate results at higher redshifts we have performed a new simulation with 288^ particles 
each of mass 2.9 x 10^ Mq with Plummer softening distance 20 kpc. In both cases the 
energy conservation, as measured by integrating the Layzer-Irvine equation, was much 
better than 1 percent. 

The comparison of power spectra in Figure 1 shows qualitative agreement between the 
second order and N-body results — in both the small dip in the spectrum at small k and 
the enhancement at high k. At early times the agreement of the two nonlinear spectra is 
excellent. This agreement extends beyond the naive regime of validity of the second order 
results. As suggested above, the dominance of the contribution from long-wave modes to 
the nonlinear enhancement at early times apparently extends the regime of validity of the 
second order results. 

At late times (0.5 ^ a ^ 1) the second order results at high k show a larger enhancement 
of the spectrum in comparison to the N-body results. There is a significant discrepancy 
in the two results even within the expected regime of validity of the second order results. 
This discrepancy, coupled with the good agreement at early times, indicates that even in 
principle the second order spectrum could not have agreed with the N-Body spectrum 
shown in Figure 1 at all times. The simple dependence of the second order spectrum on a 
given in equation (14) is incompatible with the dependence of the N-Body spectrum on a 
for the full range of k lying in the nonlinear regime. 

A part of the discrepancy at late times could arise from the dependence of the second 
order results on the upper cutoff imposed on the integrals. The cutoff dependence is 
indeed the largest at late times: for a > 0.5 reasonable variations in kc can change the 
result typically by over 10%. Another source of disagreement could be that the N-body 
simulations arc done in a finite size box, therefore they have a small-A; cutoff. Since 
the contribution from long-wave modes is positive, excluding these modes could cause 
simulations to underestimate the nonlinear enhancement of power. On comparing CDM 
simulations in boxes of sides 100 and 640 Mpc we do find this to be true, but the difference 
is very small. Thus neither of the two reasons mentioned above explain the magnitude of 
the disagreement between the second order and N-body spectra. A possible explanation is 
an inadequate suppression of the second order spectrum due to collapse on small scales, 
i.e., "previrialization" (Davis & Peebles 1977; Peebles 1990). Indeed the second order 
contribution from g > A; is negative, in qualitative agreement with such a suppression, but it 
should not be surprising if the magnitude of the suppression is significantly underestimated. 
Higher order perturbative contributions may well include some of this suppression. Our 
analytic treatment neglects small-scale pressure and vorticity, which should also suppress 
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the nonlinear enhancement of power. As we mention in Section 1, so far N-body studies 
designed to test this hypothesis have concluded that small scale effects are negligible. 
However these studies have not tested different initial spectra, and they have not examined 
the power spectrum itself with as much dynamic range as our simulations provide. 



3.3. Scaling in Time 



The nonlinear power spectrum can be used to construct statistical measures of density 
fluctuations in real space. These can then be used to study the most important consequence 
of the coupling of long-wave modes: a systematic change in the variation of characteristic 
nonlinear scales with time. We proceed to do this by first defining the rms 5p/ p averaged 
on length scale R by integrating over the power spectrum with an appropriate window 
function W: 




d^kP{a,k)W^{kR) . 



(23) 



For W we shall use three different functions: a shell in /c-space, the top-hat in real space, 
and the gaussian, given respectively by. 



W^{kR) = S^ikR - 1) , 



WTK{kR) = 



3[sm{kR) - kRcos{kR)] 



WaikR) = exp 



{kRf 
{kRf 



(24a) 
(24b) 

(24c) 



In Figure 3 we plot [47r/c"^P(a, k)]^^"^^ or ^^(a) for Wd with k = r"^, to illustrate what 
we expect for the time dependence of a characteristic nonlinear scale, denoted as i?ni(a)- 
If the spectrum evolved self-similar ly then one would expect that at all a, the onset of 
nonlinear effects occurs at a scale defined by setting 47r/c^P(a, k) — constant for some value 
of the constant of order unity. This behavior is expected for power law spectra of the form 
P{k) oc /c", and has been verified in studies of N-body simulations (Efstathiou et al. 1988). 
Even though CDM-Iike spectra are not pure power laws, the simplest assumption would be 
that they show a similar behavior. However, Figure 3 shows that at early times (small a), 
the spectrum deviates from the linear one at progressively smaller values of 4TTk^P{a, k). 
This trend is even stronger for the N-body spectra. Thus already there is a hint of a 
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systematic departure of the nonlinear scaling from the conventional expectation, due to the 
variation of the spectral index n with scale for the CDM spectrum. 

This conclusion is confirmed by using the other window functions to define the 
nonlinear scale as follows. We calculate 5ij(a) using equations (24b, c) for a range of R and 
a. We then define -Rni(a) as follows: 

dnia) = 6, for R = R^,{a) , (25) 

where 5c is a constant of order unity. Since Rn\{a) is a comoving length scale, it can be used 
to define a nonlinear scale for the mass as: M^i{a) — {47r/3)pR^i{a), where p is the critical 
density today. 

In Figure 4 we have plotted Mni(a) from a — 0.04 to 1 for the gaussian and top-hat 
filters, with 5c chosen to be 1 and 1.69 for each filter. The dependence of Mni(a) confirms 
the impression conveyed by Figure 3: nonlinear enhancement is stronger at earlier times. 
While the quantitative results depend on the choice of the window function and 6c, it is 
clear that in each of the figures the slope of the second order curve is different from the 
linear curve, and this causes the relative enhancement of Mni(a) to be larger at earlier 
times. Indeed, if the normalization of the second order curves was changed (thus shifting 
them to the right) so that at a = 1 they predicted the same nonlinear mass as the linear 
curves, then all four panels would show very similar relative enhancements at early times. 

In stating quantitative results for the time-dependence of nonlinear masses we shall 
focus on the gaussian filter with 5c — I. This choice provides the most conservative 
estimates of second order effects. At a"^ — {1 + z) — (20, 10, 5, 2), Mni{a) from the second 
order spectrum is about (180, 8, 2.5, 1.6) times (respectively) larger than the hnear case. 
Figure 4 can also be used to read off the change in the redshift of nonlinearity for the desired 
mass scale due to second order effects. (Here as in the preceding figures, the linear spectrum 
is normalized so that erg = 1 at 2; = and this fixes the normalization of the second-order 
spectrum.) For example, the mass scale 10® Mq goes nonlinear at (1 + 2;) ^ 25 as opposed 
to 19 if only the linear spectrum is used; and the mass scale 10^^ Mq at (1 -|- ;2) ~ 6 as 
opposed to 5. This change in redshift is a more meaningful indicator of the nonlinear effect, 
as the change in Mni(a) is amplified due to the steepening of the spectrum at high k. 

In Figure 4 we have also shown results from the Press-Schcchter model (Press & 
Schechter 1974, hereafter PS). The PS model is a widely used ansatz for predicting the 
distribution of bound objects of a given mass at different times (section 3.4). It relies on 
the linear growth of the power spectrum, hence it is no surprise that the shape of the PS 
curve is very similar to the linear curve. Here the PS nonlinear mass is defined as the mass 
for which a fixed fraction, 0.4, of the mass in nonlinear clumps belongs to clumps of mass 
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Mni or larger. The fraction 0.4 is chosen so that the normahzation of the PS curves is close 
to that of the other curves at a = 1 — in the upper panels it is close to the nonlinear curves 
and in the lower panels to the linear curves. 

The N-body simulations can be used to define a characteristic nonlinear mass in 
many different ways. The dashed curves in Figure 4 show nonlinear masses computed 
using the power spectrum from the simulation in the same way as for the second order 
and linear spectra above (i.e., using eq. ^5|). The results are in very good agreement 
with the second order results, as expected because of the good agreement of the second 
order and nonlinear power spectra. The relative enhancement of Mni(a) over the linear 
prediction at {1 + z) = (10, 5, 2) is (11, 2.4, 1.2). By examining all four panels it can be seen 
that, independent of the parameters used, the slopes of the curves using second-order and 
nonlinear power spectra are distinctly different from those of the linear and PS curves. The 
filled triangles use a different definition of the N-body characteristic nonlinear mass and 
will be discussed in the next section. 

Our results in Figure 4 indicate that linear scaling for M^i{a) significantly 
underestimates nonlinear enhancement at high redshift. Consequently the characteristic 
masses predicted by the PS model are much smaller than the second order and N-body 
masses for z > 4, even for the choice of parameters for which they agree at late times. This 
conclusion may appear at odds with previous tests of the PS formalism made by others. 
However, no previous tests have examined the CDM model at high redshift with as much 
dynamic range as we have. As we have emphasized, mode-coupling from long waves is 
strongest for small n; for the CDM spectrum n varies with scale and approaches —3 at high 
k. It is precisely in this limit, previously untested with high-resolution N-body simulations, 
that we find the greatest departures from linear theory and the Press-Schechter model. 



3.4. Distribution of Nonlinear Masses 

The characteristic nonlinear masses defined above do not fully characterize the 
distribution of dense clumps that form as a result of gravitational instability. A better 
comparison of theory and simulation can be made using the complete distribution of masses. 

In the N-body simulation we have identified dense clumps at mean overdensity about 
200 using the friends-of-friends (FOF) algorithm with linking distance 0.2 times the mean 
interparticle separation. (The = 288^ simulation was used at 1 + 2; = 10, while the 
A^ = 144^ simulation was used for 1 + z = 5, 2, and 1.) The distribution of nonlinear clump 
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masses is very broad, so there is no unique nonhnear mass. We have chosen to define the 
characteristic nonhnear mass Mni for this distribTition as the median clump mass defined so 
that half of the mass in clumps of at least 5 particles is contained in clumps more massive 
than M„i. The 5 particle limit corresponds to M = 1.16 x 10^^ M© and 1.45 x 10^° M© 
for the = 144^ and 288^ simulations, respectively. The resulting nonlinear masses are 
denoted by the filled triangular symbols in the upper-left panel of Figure 4; for other panels 
these points would be at the same locations as in this panel. (If the PS curves were defined 
with the same lower limit for clump masses and the same value of the mass fraction, instead 
of having no lower limit and a mass fraction 0.4, they would agree more closely with the 
N-body FOF points.) It is coincidental that this definition of N-body nonhnear mass yields 
such close agreement with the analytic predictions at ^ = 0, because the broad range of 
clump masses would allow us to vary by factors of a few. The relative variations as 
a function of redshift are more meaningful. It is clear from Figure 4 that the variation of 
these N-body masses with a departs from the linear scaling even more strongly than the 
curve computed from the N-body power spectrum. Thus, nonlinear effects on the formation 
of high-redshift objects appear to be even more significant than they are on the power 
spectrum. However, the 5 particle hmit affects the FOF characteristic mass (no lower hmit 
is imposed on the PS curve), so we should make a more detailed comparison with the mass 
distribution before reaching firm conclusions. 

The PS model makes the ansatz that the formation of bound objects is determined 
by the overdensity in the linear density field. Using the gaussian distribution of the linear 
density field, this ansatz gives the comoving number density, n{M, a) , of nonlinear objects 
of mass M in the mass interval dM at expansion factor a as (Press & Schechter 1974): 



linear rms density smoothed on the mass scale M, a{M, a), growing in proportion to a. A 
popular choice for 5c is 1.69, the value of the linear density at which a spherical top-hat 
perturbation collapses to infinite density. The PS mass distribution n{M, a) has been tested 
against N-body simulations and found by other workers to work very well. Efstathiou et al. 
(1988) tested it for scale-free simulations, and several authors have tested it for the CDM 
spectrum (e.g., Carlberg & Couchman 1989). The weaknesses of such tests — particularly, 
the finite resolution of the simulations — have been recognized by these authors, but even 
so the agreement has been surprisingly good for the range of masses and redshift probed. 
Consequently, the PS model has been widely used in predicting the number density of 
objects at high redshift, or in estimating the redshift at which a given mass scale goes 
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In this equation 5c is a free parameter which can be taken to be a constant, with the 



nonlinear. 
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Figure 5 shows the cumulative mass fraction (CMF) as a function of clump mass from 
the N-body simulation and the PS prediction. The CMF is defined by 



i.e., the fraction of mass in objects of mass M or larger. In Figure 4 we defined the PS 
nonlinear mass using the condition CMF = 0.4. For the PS prediction of Figure 5 we 
have chosen the top-hat filter with 5c = 1.69. At late times, this choice gives fairly good 
agreement for the high-mass end of the mass distribution. However, we see that at early 
times the N-body mass distribution lies systematically above the PS prediction. This is in 
qualitative agreement with the results shown in Figure 4 and supports our conclusion that 
nonlinear effects on the formation of nonlinear clumps are even stronger than they are for 
the power spectrum. 

Because coupling of long waves modifies the power spectrum and therefore the rms 
density a{M,a), the failure of Press-Schecter theory to match the N-body results exactly 
does not surprise us. As an experiment we replaced ct(M, a) in equation ( p6|) using the 
second-order and N-body power spectra instead of linear theory. The resulting CMF(M) 
falls too rapidly at large M, even after 6c is increased to compensate for the nonlinear 
enhancement of density fluctuations. If the nonlinear power spectrum is used the PS 
formula gives the wrong shape for n{M, a) because it assumes a gaussian distribution of 
densities, while the nonlinear density field has a broader distribution. We have found no 
simple modification of the PS formula that can account for the systematic departures 
evident in Figure 5. Expressing an optimistic view, we note that the PS formula is accurate 
to about a factor of 2 for the CMF over the entire range shown in Figure 5. On the other 
hand, the deviations are larger for rarer objects (smaller CMF) and the sign and magnitude 
of the deviation changes systematically with a. Therefore one should use the PS formula, 
especially at high redshift and for rare objects, only with caution after calibration by 
high-resolution N-body simulations. 



We have calculated the second order contribution to the evolution of the standard CDM 
power spectrum. We believe that our results capture the dominant nonlinear contribution 
in the weakly nonlinear regime. They are consistent with N-body results in this regime 
from z = 9 to 2; ~ 1, but show a larger enhancement of the spectrum than the N-body 
results from z ~ 1 to 2; = 0. The bulk of the second order enhancement in the growth of 




(27) 



4. 



Discussion 
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the power spectrum is provided by the mode coupling of long-wave modes, especially for 
the onset of nonlinearities at high redshift. 

By analyzing the perturbative integrals we have studied the sensitivity of nonlinear 
evolution to different parts of the spectrum, and thus have probed the dynamics of the 
mode-coupling at work. We find that on scales of interest to large-scale structure in the 
universe, the dominant contribution to the weakly nonlinear evolution of most realistic 
power spectra comes from the mode-coupling of long-wave modes. Perturbation theory is 
quite adequate for estimating this contribution since the amplitude of density fluctuations 
is small for the long-wave modes. 

An important consequence of nonlinear evolution is to change the time dependence of 
the nonlinear scale Mni(a) from linear scaling: it is found to be significantly larger at high z. 
Thus objects of a given scale go nonlinear at higher redshifts than indicated by the standard 
linear extrapolation. As discussed in Section 3, this is a consequence of the variation with 
scale of the spectral index, with n ^ — 1 on the scales of interest for large-scale structure 
and n ~ — 3 on the smallest scales. We have given quantitative estimates of this effect for 
the standard CDM spectrum for different window functions and definitions of nonlinear 
scale. For a gaussian window function and (5c = 1, which provides the most conservative 
estimates, the change in the redshift factor of nonlinearity, (1 + Zni)-, is about 20% for 
10^^ Mq objects (with linear extrapolation W Zni = 5) and increases to about 33% for 
10® Mq objects (1 + Zni = 19). We have also computed nonlinear corrections using high 
resolution N-body simulations, using the power spectrum from the simulations as well as 
directly identifying bound objects. The results are in very good agreement with the second 
order predictions, especially between z ~ 4 and 10. Quantitative comparisons are provided 
in Sections 3.3 and 3.4 and in Figure 4. 

Thus the most striking implications of second order effects arc for the formation of 
nonlinear objects at high z. Theoretical studies of, for example, the first generation of 
collapsed objects, the redshift of galaxy-formation, and reionization at high-z (see e.g., 
Couchman & Rees 1986; Efstathiou & Rees 1988; Tegmark & Silk 1993) — all require 
as an input the scale of nonlinearity as a function of z. For analytical estimates this 
is invariably obtained using linear extrapolation, as for example in the Press-Schechter 
mass distribution. We have shown (Figure 5) that the Press-Schechter theory leads to a 
systematic underestimate of the abundance of high-mass nonlinear clumps at high redshift 
in the CDM model, but have not succeeded in suggesting a simple modification that works 
better. While nonlinear coupling to long waves increases the amplitude of small-scale density 
fluctuations, it also changes the probability distribution from the gaussian distribution 
appropriate in the linear regime. 
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Most realistic cosmological spectra steepen to n ~ — 3 at the smallest scales and have 
n ^ — 1 on the largest scales of interest. This is a generic feature arising from the sluggish, 
logarithmic growth of fluctuations during the radiation dominated era, thus causing the 
scale invariant spectrum with spectral index n = 1 initially to approach n = —3 on the 
smallest scales while retaining the primeval slope on scales much larger than the size of the 
horizon at the end of the radiation dominated era. Hence for different cosmological models 
the basic features of nonlinear gravitational evolution that we have studied should hold, 
although the quantitative details would depend on the values of parameters such as as, Sc, 
Q and Hq. 

The increase in redshift of collapse relative to linear theory that we have calculated 
for CDM should also occur in all reahstic spectra provided that on the scales of interest 
n decreases sufficiently rapidly with increasing k. Our results will not apply if the dark 
matter is hot, but the qualitative implications should be the same for the evolution of 
the baryonic component in a CDM- or baryon-dominated model until dissipational effects 
become important. For spectra with a very steep slope at small scales (such as in the hot 
dark matter model), second order effects may lead to a strong nonlinear enhancement which 
would drive the spectrum to a shallower slope. 

In the near future second order power spectra from theoretical models could be related 
to the power spectrum calculated from observational surveys. Indeed the shape of the best 
fit three-dimensional power spectrum computed from results of the APM survey (Baugh 
& Efstathiou 1993) shows two characteristic features of the second order CDM spectrum: 
a relatively shallow slope at small scales and a flattening of the peak of the spectrum at 
large scales. The power spectrum computed from the CfA redshift survey (Vogeley et al. 
1992) and from the 1.2Jy IRAS redshift survey (Fisher et al. 1993) had also shown the first 
feature of a shallow slope with n just below —1 at high k, but these surveys lacked the 
depth required to determine the shape of the spectrum near the peak. It will be interesting 
to see if the extended peak of the APM spectrum is a robust feature. 
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5. Figure Captions 

Fig. 1: Linear and nonlinear power spectra at expansion factors a — 0.1,0.2,0.5, and 1, 
where a — 1 corresponds to linear as — 1. The linear spectrum is given by the dotted 
curves, the corresponding second order spectrum [P{k) — a'^Pu{k) + a'^P2{k)] by the sohd 
curves and the spectrum from high resolution N-body simulations by the dashed curve. 

The spectra are all divided by to facilitate better comparison of the nonlinear effects at 
different values of a. The triangles marked on the second order spectra indicate the point 
at which a'^P2{k) = a^Pii{k): this indicates the approximate limit of validity of the second 
order results. 

Fig. 2: Contributions to P2{k) vs. logg, where q is the magnitude of the integrated 
wavevector. The two panels are for different choices of k. P2{k) is defined in equation 
(15) and is the sum of the contributions P22{k) and 2Pi3(A;). The integrand of P22{k) is 
symmetric in q and {k — q); we have chosen to associate the contribution from such a pair 
of wavevectors with the wavevector with smaller magnitude. Other choices do not alter the 
basic trend seen here, namely, that the contribution from q < k is generally positive and 
peaked at q — k/2, while that from q > k is generally negative. Moreover, a comparison of 
the plots for the two values of k shown illustrates that at higher k the positive contribution 
from small q dominates, leading to a net enhancement of small-scale power. This is due to 
the increasing amount of power at q < k for higher k, as can be seen in Figure 3. 

Fig. 3: RMS amphtude of density fluctuations vs. scale k for several expansion factors. The 
lower curves correspond to smaller a. Sohd (dashed) curves are used for the second order 
(linear) results. The second order curves are shown only for the estimated regime of validity 
shown in Figure 1. It is clear from the results at different a that the nonlinear contribution 
becomes significant at earlier a for successively smaller values of [4:nk^P{a,k)]^^'^. 

Fig. 4: Growth of characteristic nonlinear mass with time. The mass scale Mni(a) at 
which the rms Sp/p reaches a fixed value (denoted by 5c in the figures) is plotted vs. 
the expansion factor a = -\- z). Two values of Sc are used to define Mni(a): 6c = 1 
on the left and 6c = 1.69 on the right. For each 6c the rms 6p/p is computed with a 
gaussian window function for the upper panels, and with a real space top-hat for the lower 
panels. The dotted curves show Mni(a) computed using the linear spectrum Pii{k); the 
solid curves include the second order contribution for the same normalization of the linear 
spectrum. The dot-dashed curves have been computed from the N-body power spectrum 
shown in Figure 1. The dashed curves are computed using the Press-Schechter model, 
with the characteristic nonlinear mass defined as that at which a fixed fraction, 0.4, of the 
mass in nonlinear clumps is in clumps more massive than M^i- In the top-left panel the 
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symbols labelled "N-body FOF" are obtained from the N-body simulation by using the 
friends-of-friends algorithm to identify clumps of at least 5 particles, and then to define a 
characteristic mass so that half the mass in clumps is in clumps more massive than Mni. 

Fig. 5: Cumulative mass fraction (CMF) vs. clump mass M at a = 0.1,0.2,0.5,1. The 
dashed curves represent the predictions of the Press-Schechter model, while the sohd curves 
are obtained from N-body simulations. The curves are shown at different times, with the 
higher curves representing larger values of a. The N-body curves are obtained using the 
friends-of-friends algorithm with linking parameter = 0.2. 



